classdef NumericalNMRG < handle
    properties
        gases
        beams
        conditions
        
        system
    end
    
    methods
        function obj = NumericalNMRG(varargin)
        end
        
        function obj = initializeSystem(obj, varargin)
            p = inputParser; p.KeepUnmatched = true;
            p.addParameter('cellSize',             1.5, @(x) x>0);   % cm
            p.addParameter('cellTemperature',      100, @(x) x> -273.15);  % degree Celsius
            p.addParameter('magneticField',    21.5e-6, @isnumeric); % T
            p.parse(varargin{:});
            
            obj.conditions.cellSize = p.Results.cellSize;
            obj.conditions.cellTemperature = p.Results.cellTemperature;
            obj.conditions.magneticField = p.Results.magneticField;

            obj.initializeCell(varargin{:});
            obj.initializeBeams(varargin{:});
            obj.system = CellSystem.NMRCell(obj.gases, {obj.beams.pumpBeam, obj.beams.probeBeam});
        end
        
        function obj = initializeCell(obj, varargin)
            p = inputParser; p.KeepUnmatched = true;
            p.addParameter('Rb85Abundance',     0.7217, @(x) x>=0 && x <=1); % natural abundance
            p.addParameter('Rb87Abundance',     0.2783, @(x) x>=0 && x <=1); % natural abundance
            p.addParameter('N2Pressure',           450, @(x) x>=0); % Torr
            p.addParameter('Xe129Pressure',          4, @(x) x>=0); % Torr
            p.addParameter('Xe131Pressure',         14, @(x) x>=0); % Torr
            p.parse(varargin{:});
            
            import Atom.AlkaliMetal.AlkaliMetal Atom.Other.Nitrogen Atom.Noble            
            import Atom.AlkaliMetal.Transition
            import Gas.Gas Gas.Vapor

            magB = obj.conditions.magneticField;
            rb85=AlkaliMetal('85Rb'); rb85.setMagB(magB);
            rb87=AlkaliMetal('87Rb'); rb87.setMagB(magB);
            n2 = Nitrogen();
            xe129 = Noble('129Xe');
            xe131 = Noble('131Xe');
            
            T_kelvin=273.15+obj.conditions.cellTemperature;
            g={ Vapor(rb85, 'temperature', T_kelvin, 'abundance', p.Results.Rb85Abundance), ...
                Vapor(rb87, 'temperature', T_kelvin, 'abundance', p.Results.Rb87Abundance), ...
                Gas(n2,     'temperature', T_kelvin, 'pressure', p.Results.N2Pressure*Torr2Pa), ...
                Gas(xe129,  'temperature', T_kelvin, 'pressure', p.Results.Xe129Pressure*Torr2Pa), ...
                Gas(xe131,  'temperature', T_kelvin, 'pressure', p.Results.Xe131Pressure*Torr2Pa)};
            
%             nzTest = [p.Results.Rb85Abundance, p.Results.Rb87Abundance, p.Results.N2Pressure, p.Results.Xe129Pressure, p.Results.Xe131Pressure] > eps;
%             obj.gases = g(nzTest);
            obj.gases = g;
        end
        
        
        function obj = initializeBeams(obj, varargin)
            p = inputParser; p.KeepUnmatched = true;
            p.addParameter('PumpPower',          10e-3, @(x) x>=0);           % W
            p.addParameter('PumpPolarization',  [1 1j], @(x) length(x) == 2); % circularly polarized
            p.addParameter('PumpDirection',    [0 0 1], @(x) length(x) == 3); % z direction
            p.addParameter('PumpWaist',         8.0e-3, @(x) x>=0);           % radius in [m]
            p.addParameter('PumpDetuning',         0.0, @isnumeric);          % MHz
            p.addParameter('ProbePower',         10e-3, @(x) x>=0);           % W
            p.addParameter('ProbePolarization',  [1 0], @(x) length(x) == 2); % linearly polarized
            p.addParameter('ProbeDirection',   [1 0 0], @(x) length(x) == 3); % x direction 
            p.addParameter('ProbeWaist',        1.0e-3, @(x) x>=0);           % radius in [m]
            p.addParameter('ProbeDetuning',        0.0, @isnumeric);          % MHz
            p.parse(varargin{:});
            
            import Atom.AlkaliMetal.AlkaliMetal            
            import Atom.AlkaliMetal.Transition
            import Laser.AlkaliLaserBeam
            
            rb87=AlkaliMetal('87Rb');            
            obj.beams.pumpBeam=AlkaliLaserBeam(rb87, Transition.D1, ...
                                     'name',         'Pump', ...
                                     'power',        p.Results.PumpPower, ...
                                     'polarization', p.Results.PumpPolarization, ...
                                     'direction',    p.Results.PumpDirection, ...
                                     'waist',        p.Results.PumpWaist, ...
                                     'detuning',     p.Results.PumpDetuning);

            obj.beams.probeBeam=AlkaliLaserBeam(rb87, Transition.D2, ...
                                     'name',         'Probe', ...
                                     'power',        p.Results.ProbePower, ...
                                     'polarization', p.Results.ProbePolarization, ...
                                     'direction',    p.Results.ProbeDirection, ...
                                     'waist',        p.Results.ProbeWaist, ...
                                     'detuning',     p.Results.ProbeDetuning);
        end
        
        %%
        function pumpBeam = getPumpBeam(obj)
            pumpBeam = obj.system.getCompomentByName('pump').stuff;
        end
    end
    
    methods
        result = GetCellAbsorption(obj, frequency, transition, varargin)
        result = GetBeamPropagation(obj, x, y, z, rateRelax, varargin)
    end

end

